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This  article  presents  a  model  for  the  uni-axial  compressive  response  of  uni- 
directionally  reinforced  fibrous  composite.  The  model  incorporates  the  non-linear  shear 
response  and  the  failure  strain  of  the  matrix,  accounting  for  both  aspects  within  a  non-linear 
field  equation  which  governs  the  load-deflection  process.  In  addition,  the  model  considers 
the  effects  of  two  kinds  of  geometric  imperfections,  namely,  initial  fiber  ’vavines*  and 
random  fiber  spacing.  It  is  shown  that  under  uni-axial  compression  random  fiber  spacing 
may  instigate  the  formation  of  severe  transverse  loadings  on  the  fibers,  which  suggest  the 
existence  of  a  transitional  mechanism  from  micro-buckling  to  micro-kinking. 

Computational  results  arc  presented  which  illuminate  the  effects  of  several  material 
and  geometric  factors  on  the  compressive  strength  of  composites. 


*  Department  of  Engineering  Science  and  Mechanics,  the  University  of  Tennessee  at  Knoxville. 
0  Engineering  Technology  Division,  Oak  Ridge  National  Laboratory. 


1.  INTRODUCTION 

The  compressive  behavior  of  composite  materials  has  been  studied  during  the  last 
three  decades.  Nevertheless,  many  issues  still  remain  unresolved  at  the  present  time. 
Extensive  listings  of  references  on  the  subject  can  be  found  in  several  articles,  such  as  that 
by  Shuart  [1]  and  the  review  by  Camponeschi  [2].  The  latter  author  also  lists  numerous 
factors  which  affect  the  compressive  response  and  strength  of  composites.  These  factors 
include,  among  others,  local  inhomogeneities  and  defects,  constituent  properties,  and  Fiber 
waviness  —  which  will  be  addressed  in  the  mode!  presented  in  this  work. 

This  paper  concerns  failure  by  fiber  microbuckling,  accounting  for  effects  which 
were  not  considered  elsewhere.  These  involve  a  formulation  which  incorporates  the  non¬ 
linear  shear  response  of  the  polymer  within  the  field  equation,  the  accounting  for  mtertacial 
separation  between  matrix  and  Fibers  upon  reaching  the  ultimate  shear  strain  within  the 
polymer  and,  especially,  the  incorporation  of  the  statistical  variability  of  fiber  spacings  into 
the  model.  A  significant  consequence  of  the  latter  aspect  is  that  it  induces  highly  localized 
transverse  loads  on  the  fibers,  which  could  lead  to  their  failure  by  kinking. 

Following  Rosen  [3]  and  many  other  investigators  we  employ  the  Cartesian  two 
dimensional  model  shown  in  Figure  1  (a)  where  Fiber  and  matrix  regions  are  considered  as 
lamellar  zones.  This  geometric  idealization  enables  us  to  investigate  the  effects  of  non- 
linearities  and  statistical  variations  in  Fiber  spacings  in  a  tractable  manner.  It  should  be 
noted  that  the  cylindrical-geometry  models  of  Sadowsky  et  al.  [4],  Hermann  et  al.  [5], 
Lanir  and  Fung  [6],  and  Greszczuk  [7]  were  restricted  to  linear  elastic  behavior. 

Rosen’s  model  predicts  a  compressive  failure  stress  oc  =  Gm/V,,,,  where  Gm  and 
Vm  are  the  matrix  shear  modulus  and  volume  fraction,  respectively.  This  prediction  sub¬ 
stantially  exceeds  experimental  values.  Several  modifications  to  Rosen’s  model  were 
introduced  subsequently,  in  order  to  achieve  closer  correlation  with  data.  Primarily,  these 
modifications,  such  as  the  works  of  Wang  (8],  Lin  and  Zhang  [9],  Guynn  et  al.  [10],  and 
Highsmith  et  al.  [11],  considered  non-linear  matrix  shear  response  and  initial  Fiber  wavi¬ 
ness.  The  incorporation  of  the  foregoing  material  and  geometric  factors  into  the 
microbuckling  analysis  of  Rosen  [3]  resulted  in  improved  predictions  of  compressive 
strength. 

The  analyses  of  Wang  [8]  and  Lin  and  Zhang  [9]  considered  the  non-linear  shear 
behavior  of  the  resin  but  did  not  explicitly  address  the  variation  of  the  matrix  shear  strain 
parallel  to  the  fiber  direction,  which  is  incorporated  herein.  On  the  other  hand,  it  appears 
that  the  finite  element  method  employed  by  Guynn  et  al.  [10]  and  the  iterative  numerical 
scheme  utilized  by  Highsmith  et  al.  [11]  would  become  excessively  cumbersome  in 
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handling  the  more  complicated  circumstances  of  fiber/matrix  interfacial  debonding  regions 
and  non  -uniform  fiber  spacing  considered  in  the  present  article. 


2.  FORMULATION  AND  RESULTS 

Consider  a  uni-axially  reinforced  composite  which,  following  Rosen  [3],  is  repre¬ 
sented  by  a  two-dimensional  layered  array  as  shown  in  Figure  1(a).  Let  x  and  y  denote 
Cartesian  coordinates  in  directions  parallel  and  transverse  to  the  layers,  and  designate  by  2h 
the  thickness  of  a  "fiber  layer"  centered  within  a  composite  layer  of  thickness  2c.  Conse¬ 
quently  we  have  Vf  =  h/c  and  Vm  =  (c-h)/c,  where  Vf  and  Vm  are  fiber  and  matrix  volume 
fractions,  respectively. 

Following  Rosen's  premises  [3],  we  assume  that  the  external  compressive  load  N 
is  borne  entirely  by  the  fiber  region,  which  is  modeled  as  a  Bemoulli-Euler  beam,  while  the 
matrix  responds  in  shear  only.  In  addition,  we  consider  a  micro-buckling  length  L  and  an 
initial  fiber  waviness  v&(x)  with  periodicity  of  2L.  In  the  sequel,  we  let  v&(x)  = 
50cos(jtx/L),  though  this  specific  choice  is  not  essential  to  our  method.  Finally,  in  antici¬ 
pation  of  the  circumstances  which  emerge  due  to  non-uniform  fiber  spacings,  we  denote  by 
q(x)  the  distributed  lateral  load  on  the  fiber.  See  Figure  1(b). 

Focusing  attention  on  the  case  where  all  fibers  buckle  in  phase,  the  so-called 
"shear-mode"  microbuckling  (Rosen  [3],  Garg  et  al.  [12]),  we  have  the  following  familiar 
expression  for  Y?y,  the  shear  strain  in  the  matrix: 


Y5V  “ 


—1 . 

1-Vf  dx 


(1) 


In  equation  (1)  vf  denotes  the  lateral  displacement  ( in  the  y-direction  )  of  the  fiber.  In  view 
of  the  assumption  of  Bemoulli-Euler  theory,  vf,  and  thereby  also  yjfy,  depends  only  on  x. 
Considering  non-linear  shear  response  of  the  matrix,  we  write 


T?y  =  G?F(Y«y) 


(2) 


where  G?  is  the  linear  elastic  shear  modulus  at  infinitesimal  shear  strains,  in  which  case 

F(Y?V)  *  Y  5V 

The  longitudinal  strain  in  the  fiber,  e£,  under  the  combined  effects  of  compression 
and  bending  is  given  by 

dx  2  L  dx  dx  dx 

In  equation  (3),  uf  denotes  the  fiber  displacement  in  the  direction  of  x. 
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Consequently,  the  axial  displacement  at  x=LV2  is  given  by 


"Lin 


=  A  = 


dvfox2  _  ^dvp 


dx 


dx 


/tv  4. 1 

2  EA 


(4) 


As  can  be  noted  from  equations  (3)  and  (4)  the  hypothesis  that  fibers  deform  in-phase 
implies  that  uf  and  vf  are  common  to  all  fibers  regardless  of  their  spacing.  On  the  other 
hand,  equations  (1)  and  (2)  state  that  the  support  provided  by  the  matrix  varies  with  the 
fiber  volume  fraction  Vf.  These  observations  imply  the  existence  of  lateral  loads,  q=q(x), 
which  enforce  a  common,  in-phase  deformation  of  all  fibers  in  the  case  of  non-uniform 
spacing.  To  emphasize  their  dependence  of  the  spacing  c,  we  shall  write  q=q(x,c). 

Consider  an  individual  cell  of  width  2c.  The  principle  of  virtual  work  yields 


L 


at  bet  dVf  + 


I  T*V 

Jv® 


SylPy  dV"* 


-r 


q(x,c)  8vf  dx  +  N  5A  =  0 


(5) 


Substitution  of  expressions  (1M4)  into  equation  (5)  and  employment  of  integrations  by 
parts  yield  the  following  field  equation  and  boundary  conditions  for  each  individual  cell: 

(6) 

with 


dM  =  0 

0 

11 

X 

** 

CtJ 

dx 

dx3 

< 

11 

p 

dV  =  0 

at  x  =  L 

dx2 

2 

Turning  to  the  case  of  random  fiber  spacing,  let  p(c)  denote  the  probability  density 
of  the  cell  dimension  2c.  Obviously 


p(c)  dc  =  1 


(8) 


In  the  present  circumstance  the  principle  of  virtual  work  gives 


f 


l 


-f 


p(c)  I  o{  be[  dVf  +  Tjy  by jy  dV™ -  I  q(x,c)  5vfdx  +  N  8A/dc  =  0  (9) 
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Furthermore,  in  the  absence  of  external  lateral  loads,  equilibrium  in  the  direction  of  y 
requires 


(10) 


Integration-by-parts  of  equation  (9),  upon  expressing  all  variations  in  terms  of  5(^), 
gives  the  following  field  equation  for  vf 


F(r^d^)dc+N(^+^)=0  <n) 

The  boundary  conditions  for  the  case  of  randomly  spaced  fiber  remain  the  same  as  those 
given  in  equations  (7).* 

It  is  advantageous  to  further  reduce  the  order  of  the  differential  equation  given  in 
(11)  and  express  it  in  a  non-dimensional  form  in  terms  of  the  following  non-dimensional 
parameters: 

X  =  f,  Y  =  (12) 

L  dx  L 

In  addition,  the  probability  distribution  function  p(c)  can  be  converted  to  a  probability 
distribution  function  p(Vf). 

In  view  of  expressions  (12),  the  non-dimensional  form  of  Equation  (11)  reads 

^  -  I  p(' Vf)  a2(Vf)  ( 1-Vf)  F(-*~ )  d Vf  +  X2Y  =  -  X2Y„  ( 1 3) 

dx2  jo  1-Vf 

where 

‘^-vfdbirir  ■  x2=tmr  •  Y°=-CTSinnX  <14> 

The  boundary  conditions  which  accompany  the  second  order  non-linear  differential 
equation  (13)  are 


*  In  view  of  equation  (10)  it  was  possible  to  derive  differential  equation  (11)  which  is  one  order  lower  than 
that  given  in  equation  (6).  The  lower  order  equation  (11)  enables  the  determination  of  the  lateral 
displacement  vf  to  within  a  rigid  translation,  which  is  of  no  relevance  to  the  failure  mechanisms  considered 
in  this  work.  An  additional  integration  of  expression  (11)  with  respect  to  x,  further  reduces  the  order  of  the 
differential  equation,  leading  to  a  solution  which  incorporates  an  indeterminate  rigid  body  rotation. 
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(15) 


Y(0)  =  °,  g^)-0 

In  the  case  of  uniform  fiber  spacing  equation  (13)  reduces  to 

&r- a2(l-Vf)F(-X_ )  +  X2Y  =  -  X2Y0  (16) 

dx2  I'Yf 


Note  that  for  the  linearly  elastic  case  with  uniformly  spaced  fibers  equation  (16)  takes  the 
simple  form 


with  the  solution 


&L- a2  Y  +  X2Y  =  -  X2Y0 
dx2 


Y  = - GslZ. - sjn 

X2  -  7t2  -  a2 


This  corresponds  to  the  buckling  load  predicted  by  Rosen  [3],  namely  X2  =  n2  +  a2.  Note 
that  the  above  result  assumed  an  unbounded  shear  strain  in  the  matrix  (yjy  — »  °°). 

However,  if  one  considers  a  linearly  elastic  matrix  response  followed  by  an  ideally 
plastic  deformation  atyjfy  =  yP,  then  plastic  yield  begins  at  X=l/2  and  the  onset  of  plastic 
deformation  is  found  to  occur  at 


2  (re2  +  a2)  (l-Vf)Yp 
etc  +  (l-Vf)Yp 


=  X2(yp) 


(17) 


The  above  result  agrees  with  the  value  obtained  by  Steif  [13]  for  the  slippage  initiation 
load,  beyond  which  the  matrix  no  longer  supports  the  deformed  fibers. 

Case  1:  Uniformly  Spaced  Fibers  with  Bi-linear  Shear  Modulus  of  the  Matrix 

Consider  a  bilinear  shear  stress-strain  response  of  the  matrix  material,  given  by  the 
following  expression  for  F(y?y) 


G?F(y8yH 


Gp  (Yxy~Yy)  +  G?  Y  y 


GiPYxy 

Gp(  Yxy  +  Yy  )  -  G?  yy 


if  Y  xy  >  Yy 
if  -Yy<Yxy<Yy 
if  Yxy<-Yy 


(18) 
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In  equation  (18)  GJ1  and  Gg1  are  the  slopes  of  the  stress-strain  curve  in  the  elastic  and  strain 
hardening  domains,  respectively.  It  will  be  shown  in  this  section  that  the  buckling  associ¬ 
ated  with  the  response  expressed  in  equations  (18)  can  be  handled  analytically. 

For  loads  which  correspond  to  X2  which  exceeds  X2(yy)  in  equation  (17)  the  shear 
response  of  the  matrix  will  follow  the  bi-iinear  stress-strain  relation  over  a  region 
£<X<l/2,  but  will  still  remain  linearly  elastic  within  the  central  region  0<X<^.  Obviously 
^  decreases  with  increasing  X2.  Substitution  of  expressions  (18)  into  equation  (16)  gives 

-  02  Y  +  X2Y  =  -  X2Y0  at  0<X<5 
dx2 

(19) 

-  <xl  Y  +  X2Y  -  -  X2Y0  -  P2  at  %  <  X  <  1 

dx2  ^  H  *  2 

where  cite  and  ctp  defined  according  to  (14)  with  shear  moduli  G™  and  G™,  respectively, 
and 


g2  _ _ 2i3 _ <G™  -  ^ 

P  Vf  ( 1—  Vf )  Elf 


Yy 


The  boundary  and  continuity  conditions  associated  with  equations  (19)  are 


^(1)  =  0,  Y(0)  =  0,  Y(0  =  Y(V),  =  g[<$-).  Y(^)  =  -  VmYy  (20) 

The  above  conditions  correspond,  respectively,  to  the  vanishing  of  the  moment  at  x-L/2, 
and  of  the  shear  at  x=0,  the  continuity  of  shear  and  moment  at  x=£L  and  the  requirement 
that,  by  hypothesis,  i  yjpy  i  =  yy  at  x=qL.  The  live  conditions  given  in  equations  (20) 
determine  the  four  unknowns  associated  with  the  two  second  order  differential  equations 
(19),  as  well  as  the  yet  unknown  location 

Note  that  the  solution  for  Y  determines  the  displacement  vf  to  within  arbitrary  rigid 
translations  and  rotations,  which  are  determined  from  the  requirement  of  continuity  of  vf 

and  at  x=£L,  as  well  as  vf(0)  =  0  and  ^  =  0  at  x=L/2. 
dx  dx 

The  solution  to  equations  (19)  reads: 

for  0<X<£ 


Y.(X)  =  - 


sinh  KeX 
sinh  Ke^ 


-  - sin  ^  +  (l-Vf)Yy 

X2  -  7t2  -  a| 


+  -  -£&&■ - sin  7tX  (21) 

X2  -  7t2  -  a| 
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for  ^<X<l/2  and  X2  <  a£ 


Y+(X)  =  - 


cos  kp(1-2X)/2 
cos  K^j(l  -2£)/2 


ercX2 

\2-it2~  ct2 


sin  n£,  +  (l-Vf)Yy  - 


(22a) 


+  — £s-^-2 - sin  rcX  — — ^ - 

X2  -  it2  -  ajs  X2  -  ajs 


while  for  ^<X<  1/2  and  X  >  a£ 


v  cosh  Kp(l— 2X)/2  I  pitX2  t  /i  ir  \  02  1 

Y+(X)  = - E - - —  /  SflA - sm  nE  +  (1-Vf)vv £ } 

cosh  kp( l-2^)/2  \  X2  -  tc2  _  ct^  X2  -  a g  / 


(22b) 


+  ---  -£^2 - sin  nX - 1- 

X2  -  n2  -  a£  X2  -  a£ 


In  the  above  equations  tce  =  Va|  -  X2  and  Kp  *  Vla£  -  X2I. 


Equations  (21)  and  (22)  match  all  the  conditions  (20)  except  the  continuity 

dX 


=  The  latter  condition  yields  a  characteristic  equation,  upon  differentiation  of 


equations  (21)  and  (22),  which  relates  the  position  of  ^  to  the  load  parameter  X2.  This 
characteristic  equation  must  be  solved  numerically,  with  the  physically  meaningful  solution 
corresponding  to  the  lowest  value  of  X2. 

In  our  computations  we  utilized  the  constituent  properties  reported  by  Guynn  et  al. 
(10]  for  AS4/PEEK  at  21*C.  Accordingly,  we  took  Ef  =  67  GPa,  L  =  330  pm  and  50  = 
1.65  pm  and  Vf  =  0.6.  For  purposes  of  comparison  we  also  considered  additional  values 
of  Vf  in  the  sequel.  The  non-linear  shear  stress-strain  response  was  approximated  by  a  bi¬ 
linear  relationship  withG?  =  1.3  GPa,  G™  =  0.33  GPa  and  yy  =  4.2%. 

The  resulting  stress-deflection  curve*  are  shown  in  Figure  2  for  various  values  of 
Vf.  The  symbols  "+"  on  those  curves  correspond  to  load  and  displacement  values  at  onset 
of  departure  from  linearity  in  the  shear  stress-strain  response  of  the  matrix.  Such  departure 
occurs  when  I Y  Sfy  I  =  Yy  at  X  =  1/2.  Note  that  when  Vf  =  0.9  the  composite  can  carry  com¬ 
pressive  loads  which  exceed  the  level  which  cause  departure  from  linear  matrix  response. 
However,  for  Vf  =  0.3  and  Vf  =  0.6  the  stress-deflection  curves  exhibit  the  so  called  "finite 
disturbance  buckling  behavior,"  resembling  the  buckling  of  cylindrical  shells  under  uniaxial 
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compression  or  spherical  shells  under  external  pressure  (Simitscs  f  14]).  It  is  interesting  to 
note  that  for  Vf  =  0.3  and  Vf  =  0.6  the  cusps  in  the  stress-deflection  curves,  which  corre¬ 
spond  to  maximal  load  Weis  prior  to  buckling,  occur  at  magnitudes  just  above  those  which 
cause  I  yjjy  I  =  yy  af  Y  -t/2.  It  is  obvious  that  the  theoretically  predicted  cusps  for  Vf  =  0.3 
and  0.6  cannot  realized  experimentally.  Under  load  controlled  tests  the  maximal  loads 
will  be  followed  by  total  collapse  and  under  displacement  controlled  circumstances  the 
specimen  would  snap  through  to  the  lower  load  levels  along  the  vertical  dashed  lines 
shown  in  Figure  2. 

Further  insight  into  the  compressive  response  predicted  by  the  solution  to  equations 
(19)  and  (20)  is  provided  in  Figures  3  and  4.  The  dimensionless  length  £  -  ^  )  of 

the  regions  where  the  matrix  shear  strain  I  yPy  1  exceeds  the  linear  elastic  limit  yy  is  plotted 
vs.  the  applied  compressive  stress  Oc  in  Figure  3  for  fiber  volume  fractions  Vf  =  0.3,  0.6 
and  0.9.  Note  that  Cc  increases  monotonically  with  £  for  Vf  =  0.9,  but  decreases  (after 
very  slight  initial  amplifications  )  for  Vf  =  0.3  and  0.6. 

The  variation  of  the  matrix  shear  strain  y  SPy  with  the  dimensionless  distance  X  along 
the  fiber/matrix  interfaces  is  shown  in  Figure  4  for  Vf  =  0.6.  The  four  curves  in  that  figure 
correspond  to  distinct  levels  of  non-dimensional  load  X.  The  top  curve,  with  X  =  23.10 

represents  typical  linear  elastic  results,  with  I  y  Ty  I  <  Yy  for  all  X  and  thereby  also  ^  =  0.  In 
this  case  we  obtain  a  sinusoidal  variation  of  ySPy  which  agrees  with  earlier  results  (Wang 
[8],  Lin  and  Zhang  [9]),  namely  y^y  =  A  sinrcX  with  A  =  ettX2  /  [(l-Vf)(X2-tt2-a2)].  The 
foregoing  sinusoidal  variation  persists  until  the  onset  of  inelastic  response  at  X=l/2  which 
occurs  at  X  =  Xy  =  30.79.  This  result  is  shown  by  the  dashed  line  in  Figure  4.  The  maxi¬ 
mal  value  of  the  compressive  load,  associated  with  X  =  Xm«*  =  30.81  corresponds  to  an 
inelastic  zone  of  dimensionless  length  ^  =  0.05.  In  this  case  the  variation  of  y  Jy  with  X, 
shown  by  the  dotted  line  in  Figure  4,  is  no  longer  sinusoidal.  Beyond  £  =  0.05  values  of  X 
decrease  while  A/L  increase  according  to  Figure  2.  A  typical  circumstance,  corresponding 
to  ^  =  0. 1  and  X  =  30.23,  is  shown  by  the  solid  line  in  Figure  4. 

Case  2;  Non-Uniformly  Spaced  Fibers 
Statistical  Considerations  of  Cell-Size  Distributions 

As  noted  in  the  Introduction,  non-uniformity  in  fiber  spacing  introduces  a  new 
aspect  into  the  compressive  and  buckling  behavior  of  fiber-reinforced  composites,  namely 
transverse  internal  lateral  loads  associated  with  the  common  deformation  of  the  fibers. 
Following  the  statistics  of  spatially  distributed  data  and  the  concept  of  Voronoi  cell  tessella¬ 
tion,  as  employed  to  represent  the  spatial  distribution  of  spherical  and  cylindrical  inclusions 
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(Davy  and  Guild  [15]),  we  assume  a  cumulative  distribution  function  for  the  cell  size  2c 
described  by  a  Poisson's  point  process 

P  (Oc)  =  exp  ( -2p.c  )  (23) 

In  equation  (23)  p.  is  the  frequency  of  Voronoi  cells  in  a  unit  length,  with  a  mean  cell  size 
of  p-1.  The  above  consideration  is  subject  to  the  restriction  that  fiber  regions  cannot  over¬ 
lap,  namely  c>h  ("Gibbs  hard  core  process").  Therefore,  equation  (23)  is  modified  to 
read 

P  (Oc)  -  exp  ( -2p'(c-h) )  (24) 

Since  p.*1  is  still  the  expected  value  of  the  Voronoi  cell  size,  namely. 


One  obtains 


p1  =  E(2c)  =  - 


-f 


2c  -d-  P  (Oc)  dc 
dc 


H’  = 


l-2ph 


(25) 


Equations  (23)-(25)  can  be  expressed  in  terms  of  the  fiber  volume  fraction  Vf,  as  employed 
in  equation  (13).  Let  Vf  denote  the  average  (  "nominal"  )  value  of  the  fiber  volume  fraction 
and  2c  =  p1  the  average  length  of  the  Voronoi  cells,  then  Vf  =  h/c  =  2hp.  Consequently, 
we  have 


and 


p’  = 


Vf_ 

2h  (  1-  Vf ) 


P(Oc)  =  exp 


Vf_ 

1-Vf 


Therefore,  the  cumulative  probability  that  the  fiber  volume  fraction  Vf  exceeds  a  value  Vf  is 


P(vf  >  Vf)=  1  -  P  (Oc)  =  1  -  exp 


vT1) 


The  probability  density  distribution  which  corresponds  to  equation  (26)  is 


1-  Vf '  vf  I 


(26) 


(27) 
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Computational  results  for  p(Vf)  vs.  Vf  are  shown  in  Figure  5  for  three  nominal  (average) 
values  of  Vf  ( Vf  =  0.3,  0.6,  and  0.9  ). 


The  Compressive  Response  with  Randomly  Spaced  Fibers 

The  probability  density  p(  Vf)  given  in  equation  (27)  was  incorporated  into  the  for¬ 
mulation  expressed  in  equations  (11)  and  (13)  and  employed  to  predict  the  compressive 
response  of  Gr/PEEK  (APC-2)  composite  with  Vf  =  0.6  at  a  temperature  of  T  =  21*C. 
Based  upon  the  data  of  Guynn  et  al.  [10],  the  nonlinear  shear  behavior  of  the  PEEK  resin 
was  fitted  by  a  Ramberg-Osgood  expression 


(28) 


where  GJJ1  =  1.3  GPa  as  in  the  previous  section,  A  =  94.4  MPa  and  n  =  0.12.  In  addition, 
we  took  e  =  S^/L  =  1/200  as  before  and  assumed,  somewhat  arbitrarily,  resin  failure  to 
occur  at  yjfy  =  yu  =  10%.  The  latter  assumption  was  guided  by  the  observed  tensile  failure 
at  eu  ~  4%-5%  for  PEEK  at  room  temperature  reported  by  Johnson  et  al.  [16].  The  shear 
stress-strain  response  considered  in  the  foregoing  representation  is  shown  in  Figure  6. 

The  solution  to  equation  (13),  with  Y(0)  =  0,  ^ (j)  =0*  together  with  (27)  and 

(28)  was  obtained  numerically.  Note  that  equation  (28)  was  supplemented  by  =  0  for 
tyJfy  *  >  Yu  •  To  implement  the  numerical  solution,  the  field  equation  (13)  was  expressed  by 
finite  differences  r17],  and  solved  iteratively  by  a  quasi-linearization  method. 

In  the  above  implementation,  the  probability  distribution  function  of  the  Voronoi 
cells,  p(Vf),  was  evaluated  at  100  equally  spaced,  discrete  value  of  Vf  varying  between  Vf  = 
0  and  Vf  =  1.0.  With  the  exception  of  Figure  12,  all  computations  were  performed  for 
Vf  =  0.6. 

Further  details  of  the  numerical  schemes  are  given  in  the  appendix. 

Upon  attaining  convergence  to  a  prescribed  degree  of  accuracy,  the  computational 
program  gives  the  values  of  vf,  Y,  Y  and  Y*,  as  well  as  the  shortening  of  the  column  A. 
Results  for  the  non-dimensionalized  lateral  deflection  vf/L  and  for  the  slope  Y  vs.  X  arc 
shown  in  Figures  7  and  8  for  three  values  of  non-dimensional  compressive  loads  X, 
namely  X  =  10,  20  and  26.4.  The  latter  value  corresponds  to  the  buckling  load,  since  no 
equilibrium  configuration  could  be  computed  for  X>26.4.  The  variation  of  ylfy,  the  shear 
strain  in  the  matrix  vs.  the  distance  X  at  X=26.4  is  shown  by  the  solid  line  in  Figure  9. 
This  variation  is  contrasted  with  the  variation  of  Y8V  vs-  X  for  uniformly  spaced  fibers  at 
the  same  load  level,  as  shown  by  the  dashed  line,  and  against  the  variation  of  y(fy  vs.  X  for 
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uniformly  spaced  fibers  at  X=29.5,  which  is  the  maximal  load  level  attained  in  the  uni¬ 
formly  spaced  case,  as  shown  by  the  dotted  line.  All  the  plots  in  Figure  9  correspond  to  Vf 
=  0.6  (in  the  case  of  random  spacing  Vf  =  0.6  and  the  results  are  plotted  for  the  cell  with 
Vf  =  0.6). 

Substitution  of  the  numerically  obtained  solution  for  vf  into  equation  (6)  determines 
the  lateral  load  q(x)  for  each  Voronoi  cell,  as  specified  by  its  fiber  volume  fraction  Vf. 
Results  for  q  vs.  the  non-dimensional  distance  X  =  x/L  arc  shown  in  Figure  10  for  a  typical 
"matrix  rich"  cell,  with  Vf  =  0.25,  at  load  levels  corresponding  to  X  =  10, 20  and  the  buck¬ 
ling  value  A.  =  26.4.  Similar  plots  arc  shown  in  Figure  1 1  for  a  "matrix  poor"  Voronoi  cell, 
with  Vf  =  0.95.  Note  that  sufficiently  low  levels  of  X,  i.e.  X  =  10,  yield  small  values  of 
lateral  load  q,  while  increasing  levels  of  X  raise  the  magnitude  of  q.  It  is  especially  interest¬ 
ing  to  note  the  "spikes"  in  the  plots  of  q  vs.  X.  These  localized  amplifications  occur  at 
places  where  Yxy  attains  its  ultimate  value  yu  at  some  Voronoi  cell,  with  the  sharpest  spike 
located  near  the  place  where  hy5Py  I  =  Yu  at  the  Voronoi  cell  under  consideration.  For 
instance,  the  spikes  in  q(X)  for  X  =  20  in  Figure  10  occur  at  X  =  0.15  and  X  =  0.3,  which 
are  the  locations  where  lySy  I  =  Yu  at  the  Voronoi  cells  of  fiber  volume  fractions  Vf  =  0.99 
and  Vf  =  0.98,  respectively,  at  X  =  20.  (Obviously,  the  matrix  material  in  those  cells  failed 

over  the  ranges  of  0.15  <  X  <  0.5  and  0.3  <  X  <0.5,  respectively).  On  the  other  hand,  the 
sharp  spike  at  X  =  0.25  for  X  =  26.4  in  Figure  1 1  is  associated  with  attaining  its  ulti¬ 
mate  value  Yu  within  the  very  same  Voronoi  cell  (with  Vf  =  0.95)  considered  in  that  figure, 
while  the  remaining  peaks  are  associated  with  shear  failures  in  other  cells.  Peaks  which 
occur  at  locations  X  <  0.25  are  due  to  failures  in  cells  with  values  of  Vf  >  0.95,  while 
spikes  located  at  X  >  0.25  are  due  to  failures  within  more  resin-rich  Voronoi  cells.+ 

Comparison  between  Figures  10  and  11  shows  that  resin-rich  Voronoi  cells  are 
subjected  to  relatively  lower  lateral  loads.  This  observation  is  attributable  to  the  fact  that  the 
above  mentioned  cells  sustain  shear  strains  Yxy  of  comparatively  smaller  magnitudes. 

Predicted  axial-stress  axial-strain  relations  and  compressive  strengths  under  mono- 
tonically  increasing  compressive  loads  are  illustrated  in  Figure  12  for  various  values  of  Vf. 
The  continuous  lines,  terminating  at  points  which  corresponds  to  failure,  correspond  to 
uniformly  spaced  fibers,  while  symbols  represent  computational  results  for  the  case  of 

*  It  may  seem  that  lateral  equilibrium  is  not  satisfied  for  the  individual  Voronoi  cells  since  ^qPQdX  *  0 

in  the  plots  shown  in  Figures  10  and  11.  However,  due  to  symmetry  about  X  =  0  and  X  *  0.5,  ^  q(X)  dX 
indeed  vanishes. 


12 


randomly  spaced  Voronoi  cells  with  filled  symbols  representing  failure.  The  stress-strain 
responses  shown  in  Figure  12  are  dominated  by  the  last  term  on  the  right  side  of  equation 
(4)  and  thus  remain  nearly  linear  until  failure.  Note  the  reduced  levels  of  failure  stresses 
and  strains  for  random  fiber  spacings. 

Unlike  the  circumstance  of  uniformly  spaced  fibers  with  bi-linear  shear  response  of 
the  matrix,  the  computational  scheme  for  randomly  spaced  fibers  cannot  be  extended  to 
predict  post-buckling  behavior  such  as  shown  in  Figure  2.  The  specific  values  of  the  com¬ 
puted  compressive  failure  stresses  are  listed  in  Table  1 .  That  Table  exhibits  the  effects  of 
the  nominal  volume  fraction  Vf,  the  amplitude  of  geometric  imperfection  5</L,  and  the 
presence  or  absence  of  an  ultimate  value  of  matrix  shear  strain  yu.  The  Table  also 
illuminates  the  effect  of  random  fiber  spacing. 

3.  CONCLUDING  REMARKS 

This  article  presented  a  mechanics  model  for  the  compressive  response  and  failure 
of  uni-directionally  reinforced  polymeric  composites  loaded  parallel  to  the  fiber  direction. 
The  model  accounted  for  the  non-linear  shear  response  of  the  resin,  including  its  ultimate 
shear  strain,  and  incorporated  two  kinds  of  geometric  imperfections,  namely,  initial  fiber 
waviness  and  random  fiber  spacings.  Heretofore,  the  latter  kind  of  imperfection  has  not 
been  considered  elsewhere. 

Unlike  earlier  works,  it  was  shown  herein  that  a  proper  accounting  for  the  non¬ 
linear  shear  response  of  the  matrix  yields  a  non-linear  field  equation  for  the  compressive 
behavior  of  the  composite.  In  general,  the  above  equation  could  be  solved  numerically  up 
to  failure.  Nevertheless,  in  some  special  circumstances,  it  was  possible  to  generate  a 
solution  into  the  post  buckling  range. 

Both  kinds  of  geometric  imperfections,  initial  fiber  waviness  and  random  fiber 
spacings,  were  shown  to  substantially  reduce  the  compressive  strength  of  the  composite. 
However,  random  fiber  spacings,  when  combined  with  the  foregoing  non-linear  shear 
response  of  the  matrix,  was  shown  to  introduce  imbalances  in  the  support  furnished  by  the 
matrix  against  fiber  microbuckling  —  resulting  in  highly  localized  transverse  loads  on  the 
fibers.  The  emergence  of  these  transverse  loads  could  explain  the  transition  from  failure 
through  microbuckling  to  the  more  commonly  observed  collapse  by  microkinking. 
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m 

60/L 

Uniform 

Spacing 

Random 

Spacing 

Yu  =°° 

Yu  =0.1 

Yu  =°° 

Yu  =0.1 

1360 

1360 

1381 

1292 

0.3 

0.0050 

1103 

1103 

1116 

1029 

0.0075 

941 

941 

947 

867 

2023 

2023 

2144 

1746 

0.6 

0.0050 

1541 

1541 

1583 

1234 

0.0075 

1253 

1253 

1281 

969 

4228 

4228 

4228 

2927 

0.9 

2702 

2702 

2685 

1700 

2023 

2023 

2023 

1194 

Table  1.  Comparison  of  Failure  Strength  (MPa) 
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APPENDIX;  THE  NUMERICAL  SCHEME 


The  nonlinear  second  order  differential  equation  (13)  can  be  expressed  as 

Y"  =  Q(X,Y)  (al) 

where  the  prime  denotes  derivatives  with  respect  to  X,  and 


Q(X,Y)  =  J  p(Vf)  a2(Vf)  (1-Vf)  F(^- )  dVf-  X2Y  -  X2YC 

An  error  quantity  at  i-th  iteration  step  is  defined  as 

s=  Y”(i>  -  Q(X,Y^) 

Consequently,  upon  employing  a  Taylor  series  expansion,  the  subsequent  error  quantity  is 
given  by 


<f‘iH)  =  Y"«)  (a2) 


d<j> 
dY 

in  equation  (a2) 

Y”(i+1)  -  j^r)(l)y(i+1)  =  Q(X,Y<r))  -  (^)(0Y(i)  (a3) 

Expression  (a3)  is  a  linear  ordinary  differential  equation  for  Y<i+,>  involving  the  known 
results  of  the  previous  iteration  Y<‘).  Note  that  the  derivative  of  Q  with  respect  to  Y  is 


Noting  that 


|  =  -  ^  ^  *  and  ^n-j  =  1,  we  obtain,  upon  imposing  <j> =  0 (,+1)  =  0 


f£  =  j  p(Vf)  a2(Vf)  FjyX-jdVf-X2 
Furthermore,  upon  employment  of  the  Ramberg-Osgood  model,  we  have 


F  = 


_L 


,  G?  ,  ,0-") 


Obviously,  the  boundary  conditions  in  equation  (15)  must  be  satisfied  in  every  iteration 
step. 


The  linear  differential  equation  (a3)  is  solved  by  finite  difference  scheme  as 
follows.  Divide  the  abscissa  0<X<1  into  N  equal  intervals  of  length  h=l/N.  Then  at  each 
node  X=X„=nh  the  second  order  derivative  Y"  is  expressed  as 
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Y"„  =  -^r (Yn+i  -  2Y„  +  Y„_j) 

Using  the  above  relation,  equation  (a3)  can  be  converted  to  an  algebraic  equation  of 
the  form 

*£!’  +  bf’YS*1’  +  y£|>  =  i«*»  (a4) 

Here, 


The  boundary  conditions  in  finite  difference  scheme  are  =  °  and 

The  system  of  equations  (a4)  can  be  represented  as 

A<i+t)Y(‘+l>  =  S(i+1) 

where 


(a5) 


"b<i+1)  1 

0 

yti+l)' 

/,( i+D\ 

rl 

1  b<i+,)  1 

| 

Y§+1> 

r(»+l) 

r2 

1  b<i+,) 

l 

j 

A<i+1>  = 

y(i+l)_  / 

|  S(i+,)  - 1 

. 

> 

1 

0 

2 

•r’J 

V<i+D 

\yn 

r(i+D 

rN 

Equation  (a5)  can  be  solved  by  means  of  the  LU  decomposition[17].  Accordingly,  the 
matrix  A(i+1)  is  decomposed  into  the  product  A(i+1)  =  L(i+,)U(i+,). 

Here, 
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L<i+1>  = 


0 

1  Y(ji+1)  0 

1  P(2i+1) 

1  y<*+1> 

• 

• 

i  pSTi" 

1  v(ifl) 

1  Yn-1 

0 

2  pg+"_ 

0  1 

and  P?+1)  =  bji+1) 

Pn+1)Yn+1)  =  1  (  n  =  1,2...  N-l) 

(n  =  2,3.. .N-l) 

pS+i>=bS+i)-2Yii:;) 

Denoting 

Z(i+1)  _  u<i+i)  Y<i+1)  (a6) 

equation  (a5)  is  transformed  to  L(i+1)Z(,+1)  =  S(i+1),  where  the  components  of  Z<i+1)  are 


computed  by 

•8*11 « (4w>  -  4"l”) /  PS*l)  <»  -  2,3. ..N-l ) 

The  recursive  relations  between  z^i+1),s  and  Y*,+1>,s  are  obtained  from  equation  (a6)  as 

Y<i+D_  _(i+l) 

tN  -  ZN 

Y<ni+,)  =  4i+!)  ~  Yn+1)Y<ni+11)  (n=  N-1,N-2,...,1) 

The  values  of  Y^l+1*  express  the  solution  to  equation  (13)  at  the  (i+i)th  iteration  step. 


When  X  l^n+1)  *  Ynl 2  attains  a  constant  value  within  a  prescribed  tolerance,  the  iteration 
n=l 

is  halted  and  post-processed  to  compute  deflection,  shear  strain  and  stress,  lateral  stress 
and  other  quantities. 
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Figure  3.  The  dimensionless  length,  £  =  0.5  -  £,  of  the  inelastic  zone  of  matrix  shear 

response  (|Yxyj  >  Yy  in  equation  (18))  vs.  applied  compressive  stress  for  Vf  = 
0.3,  0.6  and  0.9. 

Figure  4.  The  variation  of  the  matrix  shear  strain  yjfy  vs.  the  non-dimensionalized  distance 
X  along  the  fiber/matrix  interface  at  several  values  of  non-dimensionaiized 
applied  compressive  stress  X.  Fiber  volume  fraction  Vf  =  0.6.  Onset  of 
departure  from  linear  elastic  matrix  shear  response  at  X  =  Xy  =  30.79,  maximum 
compressive  stress  at  X  =  Xma*  =  30.81. 

Figure  5.  Distribution  of  local  fiber  volume  fraction  for  randomly  spaced  fiber  composites 
with  average  fiber  volume  fraction,  Vf,  of  0.3,  0.6  and  0.9. 

Figure  6.  Shear  constitutive  relation  of  PEEK  at  21*C  based  on  Guynn's  estimation  [10] 
with  shear  failure  strain  assumed  at  10%. 

Figure  7.  Non-dimensionalized  deflection,  vf/L,  vs.  X  for  randomly  spaced  fiber 
composite  with  Vf  =  0.6,  under  compressive  loads  corresponding  to  X  =  10,  20 

and  26.4.  Failure  shear  strain  yu  is  10%,  and  X  =  26.4  is  the  compressive 
strength  of  the  composite. 

Figure  8.  Solution  Y  of  the  governing  equation  for  randomly  spaced  fiber  composite  with 
Vf  =  0.6,  under  compressive  loads  corresponding  to  X  =  10,  20  and  26.4. 

Failure  shear  strain  yu  is  10%,  and  X  =  26.4  is  the  compressive  strength  of  the 
composite. 

Figure  9.  Comparison  of  matrix  shear  strain  for  the  Voronoi  cell  with  Vf =  0.6  in  randomly 

spaced  fiber  composite  under  its  failure  load  X  =  26.4  with  matrix  shear  strain 
for  uniformly  spaced  fiber  composite  under  the  same  load  level  and  its  own 

failure  load  X  =  29.5.  Vf  is  0.6  for  both  cases  (RS  and  US  imply  randomly  and 
uniformly  spaced  fiber  composite,  respectively). 

Figure  10.  Lateral  stress  q(X)  vs.  X  on  a  Voronoi  cell  with  Vf  =  0.25  in  randomly  spaced 
fiber  composite  with  Vf  =  0.6  at  various  levels  of  non-dimensional  compressive 

loads  X.  The  load  X  =  26.4  corresponds  to  the  failure  strength  of  the 
composite. 
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Figure  11.  Lateral  stress  q(X)  vs.  X  on  a  Voronoi  cell  with  Vf  =  0.95  in  randomly  spaced 
fib- .  composite  with  Vf  =  0.6  at  various  levels  of  non-dimensional  compressive 
loads  X.  The  load  X  =  26.4  corresponds  to  the  failure  strength  of  the 
composite. 

Figure  12.  Dimensionless  displacement  -A/L  at  X=0.5  vs.  applied  compressive  stress. 

(Solid  lines  are  for  uniformly  spaced  fiber  composite.  Symbols  are  for 
randomly  spaced  fiber  composite.  The  ends  of  lines  and  the  filled  symbols 
indicate  compressive  failure  strength  for  uniform  and  random  spacings, 
respectively.) 


Figure  1.  (a)  A  fiber  reinforced  composite  modelled  as  a  two  dimensional  lamellar  region 
consisting  of  fiber  and  matrix  plates;  (b)  a  deformed  single  cell. 
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—  A  /  L  (%) 

Figure  2.  The  scaled  compressive  displacement  A/L  at  X=0.5  vs.  applied  compressive 
stress  ac  for  various  fiber  volume  fractions  Vf  (symbol  "+"  corresponds  to  the 

circumstance  of  [y?y(X=i-)j  =  yy). 
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Figure  3.  The  dimensionless  length,  %  =  0.5  -  £,  of  the  inelastic  zone  of  matrix  shear 

response  (|YxVl  >  Yy  in  equation  (18))  vs.  applied  compressive  stress  for  Vf  = 
0.3,  0.6  and  0.9. 
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Figure  4.  The  variation  of  the  matrix  shear  strain  yty  vs.  the  non-dimensionalized  distance 
X  along  the  fiber/matrix  interface  at  several  values  of  non-dimensionalized 

applied  compressive  stress  X.  Fiber  volume  fraction  Vf  =  0.6.  Onset  of 
departure  from  linear  elastic  matrix  shear  response  at  X  =  Xy  =  30.79,  maximum 
compressive  stress  at  X  =  X™**  =  30.81. 


Figure  5.  Distribution  of  local  fiber  volume  fraction  for  randomly  spaced  fiber  composites 
with  average  fiber  volume  fraction,  Vf,  of  0.3,  0.6  and  0.9. 
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Figure  6.  Shear  constitutive  relation  of  PEEK  at  21*C  based  on  Guynn's  estimation  (10] 
with  shear  failure  strain  assumed  at  10%. 
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Figure  7.  Non-dimensionalized  deflection,  vtyL,  vs.  X  for  randomly  spaced  fiber 
composite  with  Vf  =  0.6,  under  compressive  loads  corresponding  toX=  10,  20 

and  26.4.  Failure  shear  strain  Yu  is  10%,  and  X  =  26.4  is  the  compressive 
strength  of  the  composite. 
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Figure  8.  Solution  Y  of  the  governing  equation  for  randomly  spaced  fiber  composite  with 
Vf  =  0.6,  under  compressive  loads  corresponding  to  X  =  10,  20  and  26.4. 

Failure  shear  strain  ya  is  10%,  and  X  =  26.4  is  the  compressive  strength  of  the 
composite. 


Figure  9.  Comparison  of  matrix  shear  strain  for  the  Voronoi  cell  with  Vf  =  0.6  in  randomly 
spaced  fiber  composite  under  its  failure  load  X  =  26.4  with  matrix  shear  strain 
for  uniformly  spaced  fiber  composite  under  the  same  load  level  and  its  own 

failure  load  X  =  29.5.  Vf  is  0.6  for  both  cases  (RS  and  US  imply  randomly  and 
uniformly  spaced  fiber  composite,  respectively). 
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Figure  10.  Lateral  stress  q(X)  vs.  X  on  a  Voronoi  cell  with  Vf  =  0.25  in  randomly  spaced 
fiber  composite  with  Vf  =  0.6  at  various  levels  of  non-dimensional  compressive 

loads  X.  The  load  X  =  26.4  corresponds  to  the  failure  strength  of  the 
composite. 


Figure  11.  Lateral  stress  q(X)  vs.  X  on  a  Voronoi  cell  with  Vf  =  0.95  in  randomly  spaced 
fiber  composite  with  Vf  =  0.6  at  various  levels  of  non-dimensional  compressive 

loads  X.  The  load  X  =  26.4  corresponds  to  the  failure  strength  of  the 
composite. 
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Figure  12.  Dimensionless  displacement  -A/L  at  X=0.5  vs.  applied  compressive  stress. 

(Solid  lines  are  for  uniformly  spaced  fiber  composite.  Symbols  are  for 
randomly  spaced  fiber  composite.  The  ends  of  lines  and  the  filled  symbols 
indicate  compressive  failure  strength  for  uniform  and  random  spacings, 
respectively.) 
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